Precise response functions in all-electron methods: Application to the 
optimized-effective-potential approach 
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The optimized-efFective-potential (OEP) method is a special technique to construct local Kohn- 
Sham potentials from general orbital-dependent energy functionals. In a recent publication [M. Bet- 
zinger, C. Friedrich, S. Blhgel, A. Gorling, Phys. Rev. B 83, 045105 (2011)] we showed that un- 
economically large basis sets were required to obtain a smooth local potential without spurious os- 
cillations within the full-potential linearized augmented-plane-wave method (FLAPW). This could 
be attributed to the slow convergence behavior of the density response function. In this paper, 
we derive an mcomplete-basis-set correction for the response, which consists of two terms: (1) a 
correction that is formally similar to the Pulay correction in atomic-force calculations and (2) a 
numerically more important basis response term originating from the potential dependence of the 
basis functions. The basis response term is constructed from the solutions of radial Sternheimer 
equations in the muffin-tin spheres. With these corrections the local potential converges at much 
smaller basis sets, at much fewer states, and its construction becomes numerically very stable. We 
analyze the improvements for rock-salt ScN and report results for BN, AIN, and GaN, as well as 
the perovskites CaTiOs, SrTiOs, and BaTiOs. The incomplete-basis-set correction can be applied 
to other electronic-structure methods with potential-dependent basis sets and opens the perspective 
to investigate a broad spectrum of problems in theoretical solid-state physics that involve response 
functions. 

PACS numbers: 71.15.-m, 71.15.Mb, 71.20.-b, 31.15.E- 



I. INTRODUCTION 

Density-functional theory (DFT)^ in the Kohn-Sham 
(KS) formalism^ has developed into a standard method 
for computational calculations of electronic properties 
due to its theoretical and numerical simplicity that goes 
along with the required accuracy for a large range of 
materials. In this theory, the many-body exchange and 
correlation effects that make the microscopic quantum- 
mechanical description of many-electron systems so com- 
plicated are hidden in a formally simple local potential, 
the exchange-correlation (xc) potential. Together with 
the classical electrostatic potential created by the elec- 
tronic and nuclear charges, it forms the local effective po- 
tential for a fictitious system of non-interacting electrons, 
the KS system, whose electron density, by construction, 
coincides with that of the real system. Physical quanti- 
ties, such as the total energy, inter-atomic forces, dipole 
and magnetic moments, etc., are then calculated as func- 
tionals of the electron density. 

However, the mathematical form of the xc energy func- 
tional, from which the xc potential derives, is unknown, 
and one must resort to approximations in practice. Sur- 
prisingly, already simple approximations, such as the 
local-density (LDA)^i^ and generalized gradient approx- 
imations (GGA) give reliable results for a wide range 
of materials and properties. However, as the function- 
als have been applied, over the years, to more and more 
complex materials and properties, shortcomings have be- 
come apparent: The spurious self-interaction error, in- 
herent to the LDA and GGA functionals, leads to an 



unphysical description of localized states, which as a re- 
sult appear too high in energy and give rise to erro- 
neous hybridization effects. Second, the LDA and GGA 
xc functionals do not exhibit a derivative discontinu- 
ity at integral particle numbers leading to semiconduc- 
tor band gaps that are underestimated by 40% or more 
with respect to experiment Furthermore, they show 
the wrong asymptotic behavior when exchange and cor- 
relation effects of spatially separate but interacting parts 
of a many-electron system are investigated. 

Orbital-dependent functionals form a new class of 
xc functionals jiSiii Already the formally simple exact- 
exchange (EXX) functional,—"— which treats electron 
exchange exactly but neglects correlation altogether, 
remedies the aforementioned deficiencies of the more con- 
ventional local or semi-local functionals. It has been 
shown that the EXX functional leads to KS band gaps 
that are in much better agreement with experiments^""— 
By definition, the KS band gaps do not contain the 
derivative discontinuity^i^ which indicates that the ef- 
fect of neglecting correlation is roughly of the same mag- 
nitude as the discontinuity but of different sign<22i2^ Also, 
localized d- and /-electron states appear at larger binding 
energies compared to local or semilocal functionals, which 
is a result of the exactly compensated self-interaction er- 
ror. 

The xc potential is given by the functional derivative of 
the xc energy functional with respect to the electron den- 
sity. In the case of the conventional LDA or GGA func- 
tionals, the functional derivative translates into a deriva- 
tive of a function and is evaluated straightforwardly. In 
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the case of orbital-dependent functionals, such as the 
exact-exchange (EXX) functional, the construction of the 
local xc potential is much more involved, because these 
functionals only depend implicitly on the density through 
the orbitals, and a chain-rule must be applied to evalu- 
ate the functional derivative leading to the optimized- 
effective-potential (OEP) equationj ^^i^^ 

This integral equation involves two kinds of response 
functions, one for the electron density and one for the 
KS wave functions, as well as a nonlocal term, which 
corresponds to the Hartree-Fock potential in the case of 
the EXX functional. The former response function, on 
which we will focus in the following, describes the linear 
response of the electron density with respect to changes 
in the KS effective potential and is usually calculated by 
standard perturbation theory as a sum over occupied and 
unoccupied states. 

Thus, in contrast to the Hartree-Fock method, the 
OEP approach involves the whole spectrum of unoccu- 
pied states. Physically, the unoccupied states provide 
variational freedom for the electron density to respond 
to the changes of the effective potential. In practical ab 
initio calculations, one finds that single-particle states 
up to high energies have to be taken into account, which 
calls for a basis set that is able to yield accurate wave 
functions over a very wide energy range, spanning from 
the occupied region up to high-lying unoccupied states. 
This places a much higher demand on the basis than in 
conventional LDA, GGA, or hybrid-functional calcula- 
tions, where only the electron density must be described 
accurately, while the OEP approach requires a suffi- 
ciently accurate description of the response functions, 
as well. This is a serious issue in all-electron methods 
using so-called linearized basis sets that are optimized 
for a certain energy region by construction and involve 
atom-centered numerical basis functions that are adapted 
to the local effective potential at pre-defined energies. 
Examples are the full-potential linearized augmented- 
plane-wave (FLAPW)^"— and linearized muffin-tin or- 
bital (LMT0)^"-2^ method, but also methods that rely 
on pre-calculated and tabulated basis sets constructed 
from atomic calculations)^"— A similar problem arises 
in pseudopotential approaches due to the pseudized form 
of the potential around the atomic nuclei, which yields 
accurate states only in a finite energy region. 

In a recent publication;^ we reported an implementa- 
tion of the EXX-OEP approach within the all-electron 
FLAPW methodj^"— We employed the mixed product 
basis (MPB)^-'^ to reformulate the OEP integral equa- 
tion in terms of a matrix equation that can be solved by 
standard numerical linear algebra tools. This approach 
also enables the construction of the local EXX potential 
without shape approximations. We found that the basis 
sets are not independent: the LAPW basis must be con- 
verged with respect to a given MPB. This balance condi- 
tion is a direct consequence of the sum-over-states prob- 
lem described above. Only a highly converged LAPW 
basis, which corresponds to a large number of unoccu- 



pied states, lends the density enough flexibility to react 
adequately to the changes of the potential, thus leading 
to a well converged response function. As demonstrated 
in Ref. |23, compared to conventional LDA or GGA cal- 
culations, the number of basis functions had to be in- 
creased typically by a factor of four or five, which from a 
practical point of view, cannot be a viable approach and 
effectively restricts the method to very small systems. A 
similar behavior was observed for Gaussian basis sets<^ 

In this paper, we present a numerical correction for the 
response functions, with which the balance condition is 
achieved with a considerably smaller LAPW basis lead- 
ing to much faster and numerically more stable calcula- 
tions. Furthermore, much fewer empty bands are needed 
for the construction of the response functions. The cor- 
rection relies on the observation that the LAPW basis is 
explicitly potential dependent and optimized for a given 
effective potential. Any change in this potential (other 
than a mere constant) will render the basis sub-optimal 
or even inadequate. One way to deal with this issue is 
cranking up the LAPW basis such that it can describe the 
Hilbert spaces for the unperturbed and perturbed poten- 
tials at the same time - leading to large computational 
costs as we have seen in Ref. [2^ In this paper, we pursue 
a different approach, which improves the precision of the 
calculated response functions considerably without the 
need of using larger basis sets and at the expense of only 
a small computational overhead. We derive a numerical 
correction by taking explicitly into account the changes 
of the LAPW basis induced by the perturbations. These 
changes directly follow from the potential dependence of 
the basis functions and are calculated by solving radial 
scalar-relativistic Sternheimer equations in the mufSn- 
tin spheres. Similarly, the response of the core states 
obeys fully relativistic Sternheimer equations. Addition- 
ally, we take account of the fact that the eigenfunctions 
of a Hamiltonian represented in a finite basis set are not 
exact eigenfunctions of the Hamiltonian operator, in gen- 
eral. Both corrections vanish in the limit of an infinite, 
complete basis. As they correct for different aspects of 
the incompleteness of a finite basis set, we refer to them 
as the incomplete-basis-set correction. 

Similar corrections are employed, for example, when 
shifts of the Bloch vector are considered in k • p pertur- 
bation theory^ when the magnetic moment is rotated^^ 
and when atomic positions are varied, i.e., in calculations 
of atomic forces^Siil and phonon band structures 
While in the former cases the effective potential remains 
unchanged, it does change in the latter case in principle, 
which should give rise to corresponding variations in the 
muffin-tin basis functions. Nonetheless, this effect has al- 
ways been neglected (rigid basis approximation) arguing 
that it is small, and one takes only variations into ac- 
count that are related to the spatial displacements of the 
atom-centered basis functions. The potential-dependent 
variations in the MT basis functions are thus complemen- 
tary to the incomplete-basis-set correction known from 
force calculations and potentially improve the accuracy 
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of atomic forces, too. We intend to test this conjecture 
in a future work. 

The paper is organized as foUows. Sections Ull and Hill 
give a brief introduction into the OEP formahsni and 
the FLAPW method. After a short recapitulation of the 
EXX-OEP implementation, we develop the incomplete- 
basis-set correction in detail in Section HIT] As a practical 
example we will use the EXX functional for the calcula- 
tions. However, we note that the numerical corrections 
for the response quantities are generally applicable, ir- 
respective of the employed orbital-dependent functional. 
The incomplete-basis-set correction facilitates the con- 
struction of the EXX potential considerably, in terms 
of both computational efficiency and numerical stability, 
as we will demonstrate in Sec. |V] for the case of scan- 
dium nitride. In Sec. IVII we show results for the nitrides 
BN, AIN, and GaN, as well as the perovskites CaTiOs, 
SrTiOs, and BaTiOa. Finally, we draw our conclusions 
in Sec. EIIl 



II. THEORY 

The KS formalism of DFT relies on the mapping of the 
interacting electron system onto a noninteracting system 
of KS electrons. These fictitious particles move in an ef- 
fective potential VeS (r) that is defined in such a way that 
the electron densities of the two systems coincide. The 
equation of motion for the noninteracting KS electrons 
thus reads 



(1) 



with the KS wave functions (pnk(r) and energies e„k for 
the Bloch vector k and band index n. Here and in the fol- 
lowing we restrict ourselves to the spin-unpolarized case 
and use Hartree atomic units unless stated otherwise. 
The generalization to the spin-polarized case is straight- 
forward. 

The effective potential consists of the classical elec- 
trostatic potential created by the electronic and nuclear 
charges in the system and the xc potential 



6n{r) 



(2) 



with the xc energy functional _Exc[«]- In the case of 
LDA (GGA), the energy functional is given as a simple 
function of the electron density n{r) = 2^°^'^' |(/3„k(r)P 
(and its gradient), and the functional derivative [Eq. ©] 
is evaluated straightforwardly. However, in the more 
general case of an orbital-dependent functional, which 
is only an indirect functional of the density, the ap- 
plication of the chain rule for functional derivatives is 
requested This leads, after multiplication with the 
KS single-particle response function^ 



Xs(r,r') 



Sn{r) 



occ. 
k n 



<k(r) 



(3) 



to the OEP equation in integral form 
Xs(r,r')t^xc(r')d'r' 



4EE/7:^M?^V, (4) 



k n 



where the sums over k and n run over all states. Solv- 
ing this equation yields the optimized local xc potential 
V^c (r) • We note that additional terms occur on the right- 
hand side if ii'xc exhibits an explicit dependence on the 
KS eigenvalues, as well. 

In this work, we will employ the EXX functional 



k.q n.n' 



as a practical example, whose functional derivative with 
respect to the wave functions yields the expression 



SE^ 



^:k(r')K^^r',r)dV' 



with the nonlocal exchange potential 



(6) 



III. FLAPW METHOD 

In the all-electron FLAPW method?^"— space is par- 
titioned into nonoverlapping, atom-centered muffin-tin 
(MT) spheres and the remaining interstitial region (IR). 
The tightly bound core states are completely confined 
to the spheres and are calculated by solving the fully- 
relativistic radial Dirac equation for the spherically av- 
eraged local effective potential V^Q{r), where r is mea- 
sured from the sphere center at and a is an atom 
index. 

The valence-electron wave functions are represented by 
linear combinations of basis functions 



(8) 



G 



with the reciprocal lattice vectors G. For the 
basis functions 0i(;G(r), we employ a bi-partitioned 
representation >22ii^ plane waves in the interstitial region 
with a momentum cutoff |k-|- G| < Gmax and numerical 
functions u'fp{r)Yim{r) in the MT sphere of atom a with 
the spherical harmonics Yim{j^) and a cutoff value /^ax 
for the angular-momentum quantum number I. The two 
sets of functions denoted by p = 0, 1 are matched in value 
and first radial derivative at the MT sphere boundaries 
to yield the LAPW basis functions 
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^exp[2(k + G)-r] ifreIR 
Mr) = { Y^j- ^^Lp(k,GK,(|r-R"|)yz^(r':rRa) if r e MT(a) ^^'^ 



1=0 m=-l p=0 

with the unit-cell volume Q and the matching coefficients 



ALp(k,G) = ^I'YC^k + G) cxp[^(k + G)R°] ^ J Kp(g°), j;(|k + G|5°)] , (10) 



where p = 1 ~ p, S" is the radius of the MT sphere of 
atom a, ji{x) are the spherical Bessel functions, and the 
square bracket denotes the Wronskian 

[/(»-), 9('')] /W-^;^ d^^'-^'' ■ ^ 

In the following, we restrict ourselves to the non- 
relativistic equations. The scalar-relativistic treatment 
is deferred to Appendix [X] The radial function (r) is 
the solution of the radial Schrodinger equation 



(12) 



with a pre-defined energy parameter and the radial 
Hamiltonian 

1 6>^ /(/ + 1) 

Its energy derivative ^^(r) = 9M^Q(r)/9e° is obtained 
from 



(13) 



(14) 



The energy parameters e° are typically placed in the cen- 
ter of gravity of the Z-projected density of the occupied 
states. 

The energy derivative itj\(r) provides for variational 
freedom around these energy parameters. While states 
close to the energy parameters are thus accurately de- 
scribed, the basis becomes less adequate for states that 
are energetically further away, e.g., semicore and high- 
lying unoccupied states. For these states, we may extend 
the LAPW basis with local orbitals4^"i^ These are linear 
combinations of u°Q(r), wj\(r), and solutions of Eq. ([T2|). 
Uip{r) and p > 2, with different energy parameters ei- 
ther fixed at the semicore level or at higher energies for 
the unoccupied states. The linear combinations are such 
that they vanish at the MT sphere boundary in value and 
radial derivative. Thus, the local orbitals are completely 
confined to the MT sphere and need not be matched to 
a plane wave outside. 



IV. IMPLEMENTATION 

In Ref. [13, we showed that the OEP equation [Eq. ^] 
can be cast into an algebraic matrix equation if the quan- 
tities are formulated in terms of an auxiliary basis that 



is designed to represent wave-function products. For this 
purpose we introduced the MPB^"— which is built from 
products of LAPW basis functions, giving rise to plane 
waves exp(iG • r)/^/Q in the interstitial region and MT 
functions AIlp{r)YLM{r) in the MT sphere of atom a 
with cutoff values |G| < G^^ax ^^^^ L < Lmax, respec- 
tively. (For the present purpose, the MPB functions are 
independent of k because of the periodicity of the lo- 
cal potential.) The radial parts A/£p(r) are constructed 
from the products Uip{r)uf, p, (r) with \l ~ l'\ < L < I + I' 
and also include the atomic EXX potential. We further 
form linear combinations of these functions such that 
they are continuous in value and radial derivative at the 
MT sphere boundaries as well as orthogonal to a constant 
function. (The latter is necessary to make the density re- 
sponse function Xs invertible.) For mathematical details, 
we refer the reader to Refs. I34l436i For the present work, 
we have further incorporated the boundary condition of 
zero slope for the MPB functions at the atomic nuclei. 
This is motivated by the observation that the local EXX 
potentials always show this behavior (cf. EXX potentials 
in Refs. [50l453 ). In Sec. IVII we will demonstrate that this 
constraint improves the shape of the EXX potential in 
the immediate vicinity of the atomic nuclei. 

In this way, the OEP equation [Eq. ^] for the EXX 
functional can be expressed as 



(15) 



where / is used to index the MPB functions, 

Xs,/,/ = JjM;(r)xs(r,r')Mj(r')dVdV' (16) 
is the single-particle response matrix, and 



k n 



denotes the vector of the right hand side. Inversion of 
Eq. yields the vector V^.j- The local exchange po- 
tential is then finally given by 



K(r) = ^K,/M/(r) 



(18) 
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A. Incomplete-basis-set correction 

Both the single-particle response function [Eq. ([T6|] 
and the right-hand side of the OEP equation [Eq. ([T7|) ] 
involve the derivative d(pn\i{r) / SVcs{v') , which describes 
the linear response of the wave function </5rak(r) with re- 
spect to changes of the effective potential. Eqs. (3), (16), 
and (17) show that in our formalism these changes are 
parametrized by the MPB functions {M/(r)}. We denote 
the linear response of the wave function with respect to 
M,(r) by 



(19) 



According to first-order perturbation theory, f^ki^^) 
obeys the normalization condition 



/ ^:u(r)¥'i'J,,(r)dV = 



(20) 



and the inhomogeneous differential equation 



[H-enk]^%ir) = [e^^^,, - A//(r)J ^„k(r) (21) 

with ^nli = (v«k|M/|v«k) and H ^ -iV^ + V^sir). 
Equation (pij) , the so-called Sternheimer equation^ fol- 
lows from linearizing Eq. ([1]) with respect to changes of 
the potential. Left-multiplication with the complex con- 
jugates of all other eigenstates (pn''k.{^) {n' ^ n), inte- 
gration over space, and summing over n' yield the well- 
known expression 



(1) 



{ipn'-k\Mi\ipnl< 



-ipn'kir) (22) 



and thus 



E 



<'k(i'Oynk(rO 



Vri'k(r) 



(23) 



where the sum runs over the infinite number of eigen- 
states of H. 

As the diagonalization of Eq. ([1]) in a basis representa- 
tion yields a whole spectrum of KS wave functions (p,ik(r) 
and energies e„k, comprising the occupied and a large 
number of unoccupied states, the response is straightfor- 
wardly calculated using Eq. (j23p . and one usually em- 
ploys this equation for a numerical implementation. 

However, the number of available wave functions N is 
limited in practice - it cannot exceed the number of basis 
functions A^lapw ^ so that the sum in Eq. is trun- 
cated {n' < N < A^LAPw) leading to a loss of accuracy. 
Equation then only accounts for that part of the re- 
sponse that happens to lie in the Hilbert space spanned 
by the finite number of wave functions. As a pragmatic 



solution, one can increase the number of basis functions 
and, thus, the number of eigenstates. However, as al- 
ready mentioned, we observed that Eq. ([25)) converges 
very slowly with respect to the size of the LAPW ba- 
sis. Sufficient convergence is attained only with basis 
sets that are considerably larger than the standard one 
used in GGA or LDA calculations (up to five times larger 
even for the simple example of diamond), entailing high 
computational costs rendering in part the calculations 
impractical. 

The difficulty of convergence can be overcome by two 
distinct corrections to the linear response of the wave 
function: (i) The necessity for expanding the wave func- 
tions, in practice, into a finite, incomplete basis set im- 
plies that the wave functions (^nk(r) are not pointwise 
exact solutions of Eq. ([1]). The fact that (i? — e„k)</'nk(r) 
does not vanish identically for all r in the unit cell will 
give rise to a small correction that formally resembles the 
Pulay term in atomic-force calculations, (ii) More impor- 
tantly, though, the incompleteness of the basis leads to 
a neglect of important response effects, in particular in 
the MT spheres, where the standard LAPW basis is op- 
timized for band energies close to the occupied states, 
while for higher-lying bands the basis becomes inade- 
quate. (This is generally true for linearized methods. In 
the pseudopotential approach, on the other hand, it is the 
one-particle potential itself that is constructed for the oc- 
cupied states, making it inappropriate for high energies.) 
In short, the MT functions of the standard LAPW basis 
form a poor basis for the wave-function response. This 
is not surprising since it should be easy to find a pertur- 
bation out of the many functions M/(r) that rotates the 
resulting wave function out of the Hilbert space spanned 
by the basis functions. In fact, each of the MPB func- 
tions (except for the constant function) will have this 
effect to some extent. Having said this, the question 
arises whether it is not possible to let the Hilbert space 
itself rotate in the same way as the wave function ipnk(jc) 

does so that f^ki^) remains in the co-rotating Hilbert 
space. In other words, we seek the response of the ba- 
sis functions subject to a given perturbing potential, i.e., 
<^0kG(r)/5Kff(r')- Indeed, it will turn out that this re- 
sponse is straightforwardly calculated in the MT spheres. 

We start the derivation by letting a perturbation 
SVcs{r') act on the wave function in Eq. ([S]), which for- 
mally gives 



Sfnk{r) ^ f SzG{n,k) ,/0kG(r) 
- 7T-0kG(r) + 2G(n-,k) 



(24) 

(Local orbitals and core states will be discussed later.) 
The second term arises from the fact that the basis func- 
tions 0kG(r) depend explicitly on the effective potential 
through Eqs. ([T^ and (|14p . It contains what we have 
termed basis-function response above. We will now con- 
struct this second term explicitly and then see how it 
combines with the first term. 

As the basis functions depend on the potential only in 



6 



the MT spheres, their Unear response with respect to a Linearizing Eq. ^ gives 
change of the potential is nonzero only within the spheres. 



^ Imp 



(k,G)<V(k- 



ifr e IR 
"R") ifr e MT(a), 



(25) 



where the quantities 4>'^a.ii^)^ '^ipji'^)^ ^"iip./(k, G) 
denote the linear changes of the LAPW basis function, 
the radial function, and matching coefficients, respec- 
tively, in analogy to Eq. p9| . Here, we restrict ourselves 
to MT functions Mj{r) = Miir) with angular momen- 
tum L = 0. (For simplicity, the function Mi{r) is already 
scaled with Foo(f) = l/-\/47r.) Of course, the radial func- 
tions w°p(?') can also respond to nonspherical perturba- 
tions (L ^ 0) of the potential. Then, the linear response 
consists of a superposition oi \l — . . . ,1 -\- L functions. 
We defer this more general and more complicated case 
to a later publication and note that the case L = gives 
the most important contribution. 

The functions u^p^j (r) are obtained from linearizing 
Eqs. ([T^ and which yields Sternheimer equations 

for the radial functions 



ru 



a(l). 
lOJ ^ 



Mi{r) 



(26) 



and 



'LI 



Mi{r) ™?i(r)+r<V(r). (27) 



The variation of the energy parameter is given by the 
expectation value 



1,1 



''10 



(28) 



The scalar-relativistic versions are again deferred to 
Appendix [Xj The radial inhomogeneous differential 
Eqs. and (P7|) are easily solved by integrating from 
the origin r ~ to the MT boundary r = S°' . The re- 
sulting special solutions are not uniquely defined since 
we may always add the homogeneous solution u"p(r) or 
a multiple of it. This freedom is removed by requiring 
that UiQ{r) is normalized, which leads to the additional 
conditions 



/ 



(29) 



and 



/ 



drr'ul{r)u''£}{r).m 



The Eqs. ^ and ^ together with Eq. ^ ensure that 
Uip^j{r) vanishes for constant variations of the potential. 

As an example, we show in Fig.[IJa) the radial basis re- 
sponse functions, Uoo'^/^('') ^^'^ '^*oi'/^('')i ^'^^ ^'^ atom 
of rock-salt ScN and for Z = as obtained from Eq. (pS)) 
and ((27)) together with the corresponding perturbing po- 
tential Mjir). For comparison the conventional LAPW 
s functions UQg(r) and u^lir) are presented in Fig. [TJb). 
According to Eq. the function UQQ^j''(r) is orthog- 

onal to Woo('') ^iid, as obvious from Fig. [IJa) and (b), 
it is clearly different from Uoi(?') by more than a fac- 
tor. As a consequence, it lies outside the Hilbert space 
formed by the two basis functions and will contribute to 
the incomplete-basis-set correction. A similar observa- 
tion holds for the function u^'^^j^ {r). 

The linear change of the matching coefficients 
^impC^' G) straightforwardly obtained from differenti- 
ating Eq. pO|) with respect to the potential. This finally 
gives rise to 



A 



a(i; 

Imp 



./(k,G) 



_ . .p+i [<p(^°).Ep>-4°^p>(k,G)uff,(g°)] 

^ ' K(5-),<o(^'')] ■ ^ > 

The coefficients A'^^^ j (k, G) guarantee that the resulting 

functions 0^^^ /(r), as defined in Eq. and their ra- 

dial derivatives continuously go to zero at the MT sphere 
boundaries. 

We note that the rest of the derivation applies gener- 
ally to spherical (L = 0) and nonspherical perturbations 
{L^O). Once the (fy^^l, j{r) are constructed, linear com- 
binations with the wave-function coefficients 



G 



(32) 



yield the second term of Eq. ([M)) for variations that scale 
with Mi{r) according to Eq. p^ . 

The first term of Eq. (P^ . in the following denoted 
by ip^^ j(r), lies completely in the Hilbert space spanned 

by the LAPW basis set. Accordingly, <^|,27(r) can be 
expanded in terms of the unperturbed KS wave functions 



\</5n'k|^„w)Vn'k(r) 



(33) 
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(a) 



(b) 




Figure 1: (Color Online) (a) Normalized radial basis response 
functions Woo'^'('") [('"^d) solid line] and iioi'^^('') [(blue) 
dashed line] as obtained from Eqs. 1)26^ and (|27p for angular 
momentum I = calculated in the MT sphere of the Sc atom 
in rock-salt ScN. The perturbing potential Mi{r) is shown 
as the (green) dotted line, (b) Corresponding radial LAPW 
basis functions Wqo (r) [Eq. dnj] ](red) solid line] and ugj(r-) 
[Eq. (O] [(blue) dashed line]. 



The projection coefficient (Vn'k|'^|,2 /) is obtained by ex- 
ploiting the fact that V^ikjC'^) = '^«k.7(i")+<^ik,/('') the 



solution of Eq. After left multiphcation of Eq. pij) 

with iy9*,j^(r) {n! ^ n) and integration over space one 
yields 



(e«'k - enk) ('^«'k|<^„k / + 



k|<pik,/ + ^ik./) 



^{Vn'^.\e%~MI\vr.^.), (34) 

where we additionally allow for deviations of the calcu- 
lated wave functions </5„'k(r) from the true eigenfunctions 
of the operator H by explicitly treating I?„'k(r) = {H — 
en'k)',5n'k(i") as a nonzero quantity. Using (v3„'k|<Prik) = 
and (Ai'k|<Prikj) = leads to (Vri'k|<pikj> for n ^ n' . 
The expansion coefficient for n = n' follows from the 
normalization condition Eq. (|20p . 

By adding up ^p"^^ j (r) and i^^^^ / i^) finally end up 
with the instructive result 



E 

n'<N 



S{r-v')- ^„,k(r)^:,k(r') 



ri'<Af 



(35) 



The first term contains the usual expression from first- 
order perturbation theory [Eq. ((22| [ and a correction that 
takes into account that the wave functions are not exact 
eigenfunctions of the Hamiltonian operator H due to the 
incompleteness of the basis. We call this correction the 
Pulay term in analogy to a corresponding term - the 
Pulay force - in atomic-force calculations i^°i^^'^'* As al- 
ready discussed, the first term is inaccurate because of 
the truncation of the sum {n' < N). This inaccuracy is 
corrected by the second term that arises from the ex- 
plicit variation of the basis due to a change in the effec- 
tive potential. We call this term the basis-response (BR) 
correction. In the limit of a complete basis (with an in- 
finite number of states) the Pulay and BR term would 
vanish, and the standard perturbation theory (SPT) ex- 
pression would give the exact result. The sum over states 
in the BR part can be interpreted as a double-counting 
correction that subtracts response contributions already 
contained in the first term. 

So far, we have restricted the derivation to the aug- 
mented plane waves defined in Eq. (jH]). Corresponding 
corrections can be derived for the local orbitals and the 
core states. While the derivation for the former closely 
follows the steps already presented, the core states re- 
quire some different considerations. In contrast to the 
basis functions, they fulfill the boundary condition that 
they approach zero for r — >■ oo. This boundary condi- 
tion is used to determine the energies of the core levels, 
which in contrast to the construction of the LAPW basis 
functions are not chosen as parameters but result from 
an atomic eigenvalue problem. Therefore, we use a finite- 
difference approach: we solve the atomic eigenvalue prob- 
lems for the perturbed potentials Vcs^o{r) + ^Mi{r) and 
Vc«,o{r) - ■|M/(r) with A = 0.0001, take the difference of 
the resulting core wave functions, and divide by A, which 
directly yields the linear response of the core state. As 
the fully relativistic Dirac equation is employed for the 
core states, the finite-difference approach yields the solu- 
tion of the fully relativistic Sternheimer equation. 



Finally we use Eq. ((35|) to construct the density re- 
sponse matrix [Eq. (fTB|) [ and the right-hand side of the 
OEP equation [Eq. ((T71)[: also consider Eqs. (O and P^ . 
As a result, 
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occ. unocc. 
nk n'<N 



+ 



'4E 



(Af7(p„k|<^ik,j) 



E (A//'^„k|<^n'k)(V 



(36) 



and 

*^-4EE 



occ. unocc. 



ik Ti'<JV 
',7^(1) 



{Miipn-k\(pr, 



-{Vn'w\V^''\Vnw) 



4E 

nk 



^nk ^n'k 

(' ^«k,jl-^~e"'k|yn'k) 

^nk ^n'k 



E (^ikjl^«'k>(^n'k|K^^|^„k> 



n'<N 



(37) 



are again given as a sum over an SPT, Pulay, and BR 
term. 

In its present form, the expression in Eq. ([55)) breaks 
the Hermiticity of XsjJ because the additional term is 
formally asymmetric in the indices / and J. However, 
the numerical deviation from Hermiticity is small. To 
eliminate these inaccuracies, we take the average (xs / J + 

x:,j7)/2. 



V. PERFORMANCE OF THE IBC 

We have implemented the incomplete-basis-set cor- 
rection (IBC), as described in the previous section, in 
the Fleur program package^ which is based on the 
FLAPW method. Before showing results for the nitrides 
BN, AIN, GaN, InN, and ScN as well as the perovskites 
CaTiOa, SrTiOa, and BaTiOa in the next section, we 
first analyze in detail how the convergence properties of 
the single-particle response function [Eq. (pS)) ]. the lo- 
cal exchange potential, and the band gap are improved 
by the IBC for the example of rock-salt scandium ni- 
tride. The improvements are twofold: (1) the spherical 
response function converges with much smaller LAPW 
basis sets than before; (2) for a given LAPW basis much 
fewer unoccupied states are needed for its construction. 

Unless noted otherwise, the LAPW cutoff parameters 
Gmax — 3.8 and /max = 8 were used for the calcu- 




Figure 2: (Color Online) Convergence of the trace of Xs for 
ScN as a function of the number of local orbitals wlo per I 
quantum number (0 < ? < 4) and atom. The SPT [(green) 
dashed], BR [(blue) dotted], and Pulay term [(orange) dot- 
dashed] as well as their sum l(red) solid line] are shown, re- 
spectively. We only consider spherical MPB functions, and 
the (constant) contribution of core states is neglected. 



lations of rock-salt ScN at the experimental lattice con- 
stant of 8.50 ao [tto is the Bohr radius). The Sc Is, 2s, 
and 2p states as well as the N Is state are treated as core 
states. All other states - including the 3s and 3p semi- 
core states of Sc - are treated as valence. The Brillouin 
zone is sampled with a 4x4x4 k-point mesh. 



A. Response function 

We remind the reader that the IBC consists of two 
terms, the Pulay and the BR term, which are derived as 
a correction for the expression of standard perturbation 
theory abbreviated by SPT. In principle, the core states 
are included in Eq. ([55)) as part of the sum over the occu- 
pied states. This is reminiscent of the fact that the IBC 
not only corrects for the incompleteness of the basis, but 
also comprises the exact core-state response, which will 
give a numerically important contribution, but only a 
nearly rigid shift of the convergence curves. To simplify 
the discussion, we will leave the contribution of the core 
states out until later. 

Figure |2] shows the convergence of the trace of the 
matrix Xs,iJ [Eq- ([5^ ] as well as its SPT, Pulay, and 
BR contributions for ScN as a function of the number 
of local orbitals tilo added to the LAPW basis for each 
Im channel with < / < 4 and |to| < / (in addition 
to local orbitals already used for the semicore 3s and 
3p states of scandium). The added local orbitals are 
placed at energies in the conduction band according to 
the recipe of Ref. [2^ For simplicity, we have employed 
the same number of additional local orbitals for Sc and 
N. In this case the MPB consists of 13 spherical func- 
tions, seven at the Sc atom and six at the N atom. The 
trace tr(xs) = Xsji is restricted to these functions. 



9 







5^;£__S ° ° ° ° ° ° ° 


n n a I 




ScN 


" "1.0=6 




«LO=0 


















"LO=" 










"LO=0 




"LOT U 








- »□„ ' 

LJOQooononoEODODononononononononononononon 


1 , 

"L0= (• 









50 100 150 200 250 300 350 400 



Number of unoccupied states 

Figure 3: (Color Online) Convergence of the trace of Xs for 
ScN as a function of the number of unoccupied states A'' — 
. The SPT [(green) open symbols], BR [(blue) filled 
symbols], and Pulay term [(orange) hall- filled symbols[ as well 
as their sum [(red) open symbols[ are shown, respectively. 
Circles and squares are used to distinguish the cases tilo = 
and riLO = 6. 



The slow convergence of the SPT [dashed (green) line] is 
a direct consequence of the low flexibility of the LAPW 
basis in the MT spheres with respect to changes of the 
effective potential. The dotted (blue) and dot-dashed 
(orange) line show the corresponding behavior of the BR 
and Pulay-term, respectively. As expected, both correc- 
tions become smaller as the basis set becomes more and 
more complete toward riLo = 6. The major part of the 
correction originates from the BR term, while the Pulay 
term is significant only for tt-lo = and, even there, ac- 
counts for merely 10% of the total IBC. For tilo ^ Ij 
it rapidly approaches zero. The BR term is much more 
important numerically. In fact, the dotted (blue) and 
dashed (green) lines appear to be mirror images of one 
another, showing that the BR correction compensates 
nearly exactly for what is missing in the SPT term. The 
sum of all terms produces the solid (red) line, which ap- 
pears to be constant on the scale of the diagram. From 
the inset, which shows the curve on a much finer scale, 
we see that the variations are below 0.05 %, an accuracy 
that we could never hope to achieve without the IBC. For 
this particular case, we thus do not have to employ ad- 
ditional local orbitals (i.e., tt-lo ~ 0) for the unoccupied 
states. 

Up to now, we have discussed how the IBC affects the 
convergence of the single-particle response function with 
respect to the quality of the LAPW basis. A perhaps 
more obvious convergence parameter is the number of 
states N included in the evaluation of Eq. (|36|. All 
three terms, SPT, Pulay, and BR, involve summations 
over the unoccupied states up to the maximal band in- 
dex N . Keeping in mind that the sum in the BR part 
can be interpreted as a double-counting correction, one 
can hope that the function i^^^^J ^(r) obtained nonpertur- 
batively from direct integration of the Sternheimer equa- 



tion already contains, to a certain degree, information 
about the whole infinite spectrum of unoccupied states. 
As a matter of fact, this conjecture is substantiated by 
Fig. [3J which shows the convergence of the trace of Xs 
for two different LAPW basis sets with tt-lo = (circles) 
and nLo 6 (squares) as a function of the number of 
unoccupied states N ~ Nc\/2, where Nc\ is the number of 
valence electrons per unit cell [Nd = 16 for ScN). (The 
reciprocal cutoff value of the LAPW basis has been in- 
creased to Gniax = 5.5 to generate up to 450 bands.) 
Similarly to Fig. [21 the SPT term ] (green) open symbols] 
shows a very slow convergence with respect to the num- 
ber of unoccupied states. We note that in the riLo — 
case ](green) open circles] a large part of the response is 
actually missing in the SPT term resulting in a false con- 
vergence behavior of the curve, which seems to converge, 
but toward a wrong value. As above, the BR term ](blue) 
filled symbols] is much more important than the Pulay 
term ](orange) half-filled circles] and counterbalances the 
first term almost exactly. While the Pulay term is prac- 
tically zero in the case tilo = 6, also compare Fig. [51 
it is significant in the case tilo = and yields a small 
but important contribution to the sum of all terms. In 
fact, this sum nearly follows the same curve ](red) open 
symbols, circles and squares correspond to tilo = and 
?iLO = 6, respectively] in the two cases, showing again 
that with the IBC we can restrict ourselves to the con- 
ventional LAPW basis, i.e., tilo = 0- The total sum con- 
verges extremely fast, thanks to the BR term, in which 
the infinite spectrum of eigenstates is already incorpo- 
rated to a large extent. Without showing further results 
we note that the IBC affects the convergence of the right 
hand side ]Eq. ([57])] in a similarly beneficial way. 

So far, we have not discussed the contribution of the 
core states to the density response Xs- For the core state 
response we solve a fully relativistic Sternheimer equation 
by a finite-difference approach as discussed in Sec. IIV Al 
The resulting solution embodies the full infinite spectrum 
of unoccupied states by construction and can, therefore, 
be considered to represent already the exact core-state re- 
sponse (up to numerical errors connected with the finite- 
difference approach, which can be made arbitrarily small, 
though). Therefore, we can set to the number of occu- 
pied states in Eq. (|36p . The first term is then zero, and 
only the BR term remains. 

The contribution of the core states to the response 
function is numerically important. In the case of ScN 
it is more than four times larger than the valence con- 
tribution. However, it should be noted that the effective 
quantity for the construction of the local exchange po- 
tential is not the response function itself, but its inverse 
]see Eq. Therefore, large eigenvalues of Xs become 

comparatively unimportant in x^^- This is confirmed 
by the following observation. We find that the SPT ex- 
pression alone is incapable of describing the core-state 
response. Even with tilo = 6 and N = 450, only about 
20% of the contribution of the core states is accounted 
for, which manifests the hardly surprising fact that the 
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LAPW valence basis is unsuitable to describe changes in 
the core states. However, this shortcoming affects the 
resulting exchange potential and KS transition energies 
only slightly, as we will see below. 



B. Local exchange potential 

Figure [4] shows the local exact exchange potential be- 
tween neighboring scandium and nitrogen atoms for three 
calculations, all of them without any local orbitals for the 
unoccupied states. For these calculations an MPB with 
^max = 2.8 and Lmax = 4 is employed, and the IBC 
is applied to both core and valence electrons. Without 
the IBC (green dashed line) we obtain an unphysical, 
strongly varying potential that even tends to an unrea- 
sonable positive value at the position of the nitrogen nu- 
cleus. In Ref. [iOwe showed that the spurious oscillations 
can be avoided by augmenting the LAPW basis with lo- 
cal orbitals leading to a smooth and physical potential, 
at the expense of a very costly calculation. As shown in 
Fig. m the same smooth shape of the local exact exchange 
potential is realized by employing the IBC with a consid- 
erably smaller computational overhead. We emphasize 
that no local orbitals are used in this calculation. 

When one watches closely, one still sees very slight 
anomalies of the dotted curve at the atomic nuclei, here 
more pronounced at the nitrogen nucleus (also cf. the 
insets). These result from the fact that the radial MT 
potential enters with a factor into the equations. The 
region close to the nuclei has only little weight and is 
therefore difficult to converge. On the one hand, for small 
r the total effective potential is dominated by the poten- 
tial of the nucleus —Z/r as well as the centrifugal kinetic 
energy barrier l{l + l)/(2r^) such that slight inaccura- 
cies at the nuclei prove to be irrelevant in the calculation 
of the electronic structure. On the other hand, we can 
find a simple remedy by an additional constraint for the 
spherical functions (L = 0) of the MPB, which has the 
additional benefit of reducing their number by one for 
each atom in the unit cell. We require the gradient of 
these functions to vanish at r = 0. As already mentioned 
in Sec. IIIIl this behavior has been observed for the local 
exact exchange potential in previous publications i^"— 
Thus, the constraint does not induce errors. With such 
a modified MPB the anomalies of T4(r) at the nuclei dis- 
appear, and we obtain the red solid line shown in Fig. 2] 
We also note that the numerical stability benefits from 
this modification. 

As already pointed out, the IBC is presently only ap- 
plied to spherical variations of the potential. To converge 
the nonspherical contributions properly we still need a 
few local orbitals. Figure [SUa) shows this effect, which is 
strongest close to the atomic nucleus of Sc. Using a single 
set of local orbitals, i.e., riLO = Ij gives a slight correc- 
tion of the potential, which gives rise to changes in single- 
particle transition energies in the order of 0.10 — 0.15 cV. 
(We have considered the gap transitions F — ^ F, F — X, 
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Figure 4: (Color Online) Local EXX potential for ScN on a 
line connecting the neighboring Sc and N atoms. We have 
used a conventional LAPW basis without local orbitals for 
the unoccupied states. The (blue) dotted and (green) dashed 
curves correspond to calculations with and without the IBC. 
For the (red) solid curve we have employed an additional 
constraint for the MPB (see text). Eq. ([4]) defines the po- 
tential only up to a constant. Here, we use the convention 
JV4r) dV = 0. 

and F — >■ L.) For riLO = 2 changes are less than 0.03 eV. 
It is important to note that the calculations always con- 
verge to the same result, irrespective of whether or not 
the IBC is used. We expect that once the IBC is ex- 
tended to the nonspherical MPB functions, no extra local 
orbitals are required anymore. 

As already discussed, taking into account the exact 
core response yields a numerically important contribu- 
tion to the single-particle response function. To demon- 
strate the effect on the local exchange potential, we 
switch the IBC for the core states off. Figure [Sjb) shows 
that this has only a comparatively small effect on the 
shape of the resulting potential. At first sight, the in- 
curred changes in the potential are more pronounced than 
in Fig.[5ja), but, in contrast to Fig.[5^a), they mostly af- 
fect the region close to the nucleus, where the effective 
potential is dominated by the potential of the nucleus 
and the angular kinetic energy. In fact, the single-particle 
transition energies are influenced only little: they change 
by less than 0.04 eV. However, it should be noted that 
the exact core response is evaluated at virtually no extra 
cost. 



C. KS band gap 

After discussing the effect of the IBC on the ingredi- 
ents of the OFF integral equation and the local exact 
exchange potential, we now turn to the KS band gap. As 
we have presently derived the IBC only for the spherical 
MPB functions, we employ the cutoff Lmax = for this 
test calculation in order to highlight the improvements 
resulting from the IBC. This will produce a potential 
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Figure 5: (Color Online) (a) Comparison of I4(r) close to 
the atomic nucleus of Sc for calculations with [(green) dashed 
line, riLO = 1] and without [(red) solid line, tilo = 0[ local 
orbitals. (b) V^(r) with [(red) solid line[ and without [(green) 
dashed line[ the core response of the IBC. 
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Figure 6: (Color Online) Convergence of the direct gap of ScN 
with respect to the number of local orbitals (riLo) added to 
the LAPW basis for one-shot EXX calculations with [(red) 
solid line[ and without [(green) dashed line)[ the IBC. The 
computation time is shown on the right scale. 



Vji{r) that is purely spherical in the mufEn-tin spheres. 
We note again that a generalization to nonspherical func- 
tions (L > 0) is possible, but not yet implemented. Fur- 
thermore, we restrict the interstitial exchange potential 
to a constant for simplicity. 

Figure [5] shows the direct (F — > F) gap of ScN for 
this specific numerical setup as a function of the num- 
ber of local orbitals per I channel added to the LAPW 
basis. The case tilo = 6 corresponds to 300 additional 
basis functions. It is known that small inaccuracies in 
the band gap can occur due to the linearization error of 
the LAPW basis (in the case tilo = 0). To eliminate 
this pure basis-set effect, we have taken two measures: 
(1) we have performed only a single EXX-OEP iteration 
(starting from a PBli ground state), and (2) the final di- 
agonalization of the KS Hamiltonian was performed with 
the most accurate basis set (riLO = 6). The value riLo 
shown on the abscissa thus corresponds to the basis used 
in solving the OEP equation, and variations in the band 
gap can be attributed exclusively to the precision of the 
local exchange potential without additional basis-set ef- 
fects. We note that, in spite of the very small MPB used 
here and in spite of performing only one iteration, the 
resulting band gap is surprisingly close to the fully con- 
verged one (see Tab.|I|. 

Without the IBC [(green) dashed curve] four sets of 
local orbitals (riLo = 4) are necessary to obtain a direct 
gap with an accuracy of 0.05 cV (we note that, judging 
from the form of the curve, the accuracy at riLo = 3 
seems to be due to a fortuitous cancellation of errors). 
On the other hand, the calculation with the IBC yields a 
reliable gap with an accuracy of 0.007 eV already without 
any local orbitals for the unoccupied states (tt-lo = 0). 
Both curves converge to the same band-gap value, while 
the convergence of the IBC values is hardly visible on the 
scale of the diagram. We have indicated the computation 



time on the right scale. The computational overhead of 
the IBC calculations is due to the additional evaluation 
of the Pulay and BR term, in particular for the matrix 
elements {ipni<i\V^^\(fn-k.i) in Eq. ((37|) . It is difficult to 
compare the efficiency of the calculations, since the least 
accurate calculation with the IBC is still more accurate 
than the most accurate one without. If we take an ac- 
curacy of 0.05 eV as a criterion, one would deduce an 
acceleration of the code by a factor of four. Using less 
unoccupied states (cf. Fig. could further reduce the 
computation time. 



VI. RESULTS AND DISCUSSION 

In Table |T] we report KS transition energies, i.e., 
KS eigenvalue differences, for the III-V nitrides in the 
zincblende structure and for rock-salt ScN, calculated 
with the EXX and EXXc functionals. The latter con- 
tains additionally the LDA correlation functional in the 
parametrization of Perdew and Zunger.— While the 
ground-state crystal structure of the III-V nitrides is 
wurtzite, they can be synthesized in the metastable 
zincblende structure by epitaxial growth techniques. All 
calculations were performed at the experimental lattice 
constants (ScN, 8.50 ao; BN, 6.84 aq; AIN, 8.26 aq; GaN 
8.50ao; InN 9.41 ao) with an 8x8x8 k-point sampling 
and with local orbitals for the complete semicore shell 
(e.g., 2s and 2p states of Al). The numerical cutoff pa- 
rameters Gmax, GJnaxi ^max, and L^ax as wcll as the 
number of local orbitals are determined such that the 
transition energies are converged to within lOmeV. We 
compare our results with plane- wave pseudopotential cal- 
culations and experimental values from the literature. 

The EXX and EXXc functionals give KS transition 
energies that are much closer to the experimental value 
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than LDA. InN and ScN, which are metalhc in LDA, are 
correctly predicted to be semiconductors. The inclusion 
of the LDA correlation functional in EXXc increases the 
values by about 0.01 — 0.2 eV, but does not yield a sys- 
tematic improvement when compared to experiment. For 
AIN and ScN our EXXc values agree well with those from 
the pseudopotential studies. The differences are mostly 
of the same order as the differences in the LDA transition 
energies. 

In the case of GaN and InN the situation is more com- 
plex because of the semicore d states. In pseudopotential 
calculations they can be treated approximately as atomic 
core states in the pseudopotential or, at a considerably 
larger computational expense, as valence electrons within 
the plane-wave basis. Therefore, Table U lists two values 
for each system and functional from calculations where 
the d semicore electrons are treated as valence ("with d") 
or not ("no d"). The former values, which are systemat- 
ically smaller than the latter, must be considered to be 
more accurate. Of course, we do not have to make such 
a distinction in our calculations, because the FLAPW 
method treats valence and core electrons, down to the 
Is state, on an equal footing. Indeed, our LDA tran- 
sition energies compare better with the pseudopotential 
calculations "with d" . However, we find a much larger 
discrepancy for the EXXc functional. In all cases the all- 
electron values fall in-between the two pseudopotential 
results, often being equally far away from either value. 
Engel and Schmidt pointed out for the case of transi- 
tion metal monoxides that the complete semicore shell, 
comprising the 3s, 3p, and 3o? states, must be treated as 
valence to obtain accurate EXX values within the pseu- 
dopotential approach. It is to be expected that this also 
holds for GaN and InN, which would explain the discrep- 
ancies observed for the pseudopotential calculations with 
respect to the all-electron results. 

Table [IT] shows the energetic positions of the d levels 
in GaN and InN with respect to the Fermi level, which 
is fixed at the valence band edge. It is well known that 
the LDA underbinds the 3d states, which is usually at- 
tributed to the self-interaction error of LDA. However, 
in spite of the fact that the EXX(c) is free of this error, 
the position of the d bands is hardly improved. While 
in GaAs the d levels are lowered by about 1 eV, they re- 
main nearly at the same energy in InN. In comparison 
to the experimental values, 17.7 eV (GaN) and 14.9 eV 
(InN; measured for wurtzite structure), this is hardly an 
improvement over LDA. On the other hand, we find that 
the HE method yields rf-band positions at 21.28 cV (GaN) 
and 18.22 cV (InN), respectively, i.e.. even below the ex- 
periment. As the HE and EXX-OEP methods employ 
the same energy functional, the question arises why the 
d-band positions appear at different energies. In short, 
what is the meaning of the EXX-OEP single-particle en- 
ergies? According to Koopmans' theorem,— the single- 
particle energies in the HE approach can be understood 
as total-energy differences between an excited state and 
the ground state, neglecting orbital relaxation effects. 



For the d states, AE = E^'^ - = {d\h^^\d) = ef^, 
where the two total energies correspond to the many- 
body states with and without a hole in a d level, e^^ is 
the HE single-particle energy of that d level, and h^^ 
is the HE single-particle Hamiltonian. The OEP ap- 
proach was originally intended as a procedure to sim- 
plify the solution of the HE equations by introducing a 
local potential that minimizes the HE total energy. In 
fact, the EXX-OEP and HE ground-state total energies 
are nearly identicaL ^"'^^ This indicates that the slight 
difference in shape of the single-particle wave functions 
only has a negligible effect on the ground-state total en- 
ergy E^ , a statement that should also hold for the hole- 
state energy E^~^. Using Koopman's theorem again 
with the EXX-OEP wave functions instead of the HE 
ones, we should obtain practically the same AE. How- 
ever, as seen in Table [Hi the HE eigenvalue ejf^ differs 
from the corresponding EXX-OEP eigenvalue e^^^ by 
A = P - ef^^ = {d\V^^ -V^\d). As a consequence, 
the underbinding of the d bands in LDA - and also in 
EXX(c) - cannot be solely ascribed to the self-interaction 
error contrary to common belief. Part of this underbind- 
ing is due to the forced locality of the effective potential, 
which gives rise to a nonzero energy shift that is very 
similar in spirit to the derivative discontinuity2i2i of the 
xc functional. Even with the exact functional we cannot 
expect that the d-band positions lie at the experimental 
positions in the same way as we cannot expect the KS 
band gap to equal the experimental gap. 

Finally, we report results for the electronic structure of 
the perovskite transition-metal oxides CaTiOa, SrTiOs 
and BaTiOs. We have calculated the materials in the 
ideal cubic structure at experimental lattice constants 
(CaTiOa, 7.35 ao; SrTiOg, 7.46 ao; BaTiOg, 7.60 ao) and 
with the functional LDA, EXX, and EXXc. The Ti 3s 
and ip states as well as the first subshell of the cations 
(Ca. 3s 3p; Sr, 4s 4p; Ba 5s hp) have been treated with 
local orbitals. The electronic structures of the three sys- 
tems are very similar. They are semiconductors with an 
indirect band gap between the R and the F point and 
a direct gap at the F point. We compare the LDA and 
EXX band structures and densities of states of CaTiOa in 
Fig. El The EXX functional shifts the Ti 3d conduction 
bands away from the O 2p valence bands and thus opens 
the band gap, while the band dispersions remain nearly 
unaffected. The valence band width is slightly decreased, 
though. 

KS transition energies are reported in Table IIIII The 
LDA underestimates the energies by roughly 50%. Al- 
though the EXX(c) functionals yield values that are 
closer to experiment, we observe a pronounced overesti- 
niation by about 25%. This might be due to two reasons: 
(1) The fortuitous error cancellation arising from the ne- 
glect of correlation on the one hand and of the derivative 
discontinuity on the other hand does not work for these 
transition-metal oxides. (2) The ideal cubic structure 
used in our study does not correspond to the experimen- 
tally measured systems. It is known that these perovskite 
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Table I: KS transition energies (in eV) obtained with the LDA, EXX, and EXXc potentials for the Ill-V nitrides and ScN. An 
8x8x8 k-point sampling has been employed for all calculations. For comparison, plane-wave pseudopotential results and 
experimental values from the literature are given. 
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Table II: Energetic positions of the d levels in GaN and InN 
in eV with respect to the Fermi energy and averaged over the 
Brillouin zone and the d bands. 
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materials undergo a number of structural phase transi- 
tions. For example, SrTiOs is cubic at room temperature 
but tetragonal at low temperatures. BaTiOs, on the con- 
trary, is cubic at elevated temperatures and tetragonal at 
room temperature. 



VII. CONCLUSIONS 



We have described an efficient way to calculate pre- 
cise all-electron response functions within the FLAPW 
method. The key is the development of an incomplete- 
basis-set correction (IBC), which was derived by describ- 
ing the response of the LAPW basis functions to changes 
of the effective potential explicitly by means of radial 
Sternheimer equations. In this way, the IBC incorpo- 
rates important response contributions that lie outside 
the Hilbert space formed by the LAPW basis set. The 
resulting formula for the response function consists of 
three terms: the conventional sum-over-states expression 
of standard perturbation theory (SPT term) and two ad- 
ditional terms, the basis-response and the Pulay term, 
whose mathematical expression resembles that of the Pu- 
lay force of atomic-force calculations. The basis-response 
term is more important numerically than the Pulay term. 
Both vanish in the limit of a complete, infinite basis. To- 
gether they constitute the IBC, which also yields an ex- 
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Figure 7: (Color Online) Comparison of the LDA and EXX band structures and densities of states (DOS) for CaTiOa. The 
total and partial DOS for Ca 3d, Ti 3d, and O 2p are shown as (black) solid, (blue) dot-dashed, (green) dashed, and (red) 
dotted lines, respectively. 



Table III: Direct and indirect KS band gaps (in eV) for cubic 
CaTiOg, SrTiOs, and BaTiOa compared with experimental 
values. A 6x6x6 k-point sampling was employed. 
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plicit and, in principle, exact treatment of the response 
of the core states that avoids the sum-over-states expres- 
sion. The total correction is not small. It can be as large 
as the SPT term. 

As a practical example, we have employed the IBC to 
the EXX-OEP method, an approach that allows the con- 
struction of a local exchange potential from the orbital- 
dependent EXX functional. It involves two response 
quantities: the density response function and a response 



function for the single-particle states. We have imple- 
mented the IBC for both and demonstrated explicitly for 
the case of rock-salt scandium nitride that it improves 
the convergence with respect to the LAPW basis and 
the number of unoccupied states considerably. While 
without the correction the solution of the OEP equation 
requires a highly converged LAPW basis with a large 
number of local orbitalSf2£ no extra local orbitals are 
needed when we use the IBC. A similar statement can be 
made about the number of unoccupied states in the sum- 
over-states expression related to perturbation theory. As 
the numerical solution of the Sternheimer equation al- 
ready incorporates the infinite single-particle spectrum 
to a certain degree, the response function converges at 
much fewer bands. 

With this new scheme we have performed all-electron 
EXX(c)-OEP calculations for the III-V nitrides in the 
zincblende structure as well as rock-salt ScN. Our results 
agree favorably with plane-wave pseudopotential EXXc 
calculations from the literature. However, larger devia- 
tions are observed for GaN and InN, which we attributed 
to the neglect of semicore states (3s, 3p of Ga and 4s, 
Ap of In) in the pseudopotential calculations. This is 
in accordance with a recent publication by Engel and 
Schmid.'S Despite the fact that the EXX functional is 
self-interaction-free, it does not necessarily improve the 
position of the d states with respect to LDA, as shown for 
the semicore d states of GaN and InN. We have explained 
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this observation by the forced locality of the exchange po- 
tential, which has the effect that the KS eigenvalues of 
the occupied states cannot be associated with ionization 
energies via Koopmans theorem. In order to invoke the 
latter, extra terms similiar in spirit to the derivative dis- 
continuity of the xc potential have to be considered. Con- 
sequently, even with the exact xc potential the d states 
cannot be expected to appear at the experimental posi- 
tions. 

Furthermore, we discussed the EXX(c) electronic 
structure of the three cubic perovskites CaTiOs, SrTiOs, 
and BaTiOs in comparison with LDA. The EXX func- 
tional opens the band gap in all materials but leaves the 
dispersion of the bands nearly unaffected, except most 
notably for a slight reduction of the valence band width. 
While LDA underestimates the band gaps by about 50%, 
the EXX(c) potential overcompensates this error and 
leads to an overestimation of 25%, which we attribute 
to an incomplete error cancellation between the neglect 
of correlation and the neglect of the xc derivative dis- 
continuity. On the other hand, differences in the crystal 
structure between the theoretical setups and experimen- 
tal realizations might also play a role. 

We note that the IBC is a general approach, which 
is not restricted to the FLAPW method. It applies to 
all electronic structure methods with a basis set opti- 
mized for the effective potential, including those which 
are atom-centered and based on pre-calculated and tab- 
ulated basis sets. Typical electronic structure methods 
of this type include LMTO^-^ DM0L3 FHIAIMS^ 
OPENMX^ and the SIESTA code^i 

Furthermore, with suitable generalizations the IBC can 
be applied to a broad spectrum of response functions in 
solid-state physics, e.g., response quantities arising from 
the displacement of external potentials (e.g. phonons, 
elastic constants, stress tensor) or from external fields 
(e.g. g-tensor, chemical shift and nuclear magnetic res- 
onance, dielectric response, infrared and Raman inten- 
sities, magneto-elastic and magneto-electric tensors), or 
from more general perturbations of the Hamiltonian (e.g. 
Born effective charges, polarizability), to name a few. 

Another method that may profit from the IBC is the 
GW approximation for the electronic self-energy, which 
involves the calculation of a dynamical density response 
function and the single-particle Green function. It is well 
known, and Zinc oxide is a prominent example of it^iS 
that similarly to the KS response function in EXX-OEP, 
GW calculations usually converge badly with respect to 
the basis set and the number of unoccupied states. 
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Appendix A: Scalar-relativistic equations 

In the scalar-relativistic approximation^ for the Dirac 
equation the large and small radial component, r~^p(r) 



and (cr) ^q{r), of an electron in a spherical potential 
Vcfffi{r) at energy e obey the equations of motion 



= 2m,c\q{r) + -p{r) 
dr r 



dq{r) 
dr 



1 



q{r) + wp{r) 



(Al) 
(A2) 



with jTiici = 1 + [e — V"off,o('')]/(2c^) (c denotes the speed 
of light) and w = [l{l + l)]/[2mrcir^] -|- VeSfi{r) - e. For a 
given spherical perturbation that scales with Af (r) (index 
/ omitted to simplify the notation) , the linear changes of 
p(r) and q{r) are given by the solutions of the differential 
equations 



= 2m,eig'(r) + -p\r) 
dr r 



1 



dq'jr) 
dr 



[e' - M{r)]qir) 



-[1 



q'{r) + wp'(r) 
/(/ + 1) 



(A3) 



(A4) 



t - M{r)]p{r) 



By direct differentiation we find the differential equations 
for the energy derivatives p{r) ~ dp{r)/de and q{r) = 
dq{r)/de 



'^^^''^ - 2m,,ig(r) + ip(r) + lg(r) (A5) 



dr 
dqir) 
dr 



--q(r) + wpir) 
r 

+ 1) , , , 



(A6) 



as well as for the linear changes of p{r) and q{r) 



dp'(r) ^ ... . 1 .,, , 

dr r 
1 



(A7) 



+ -{[e'-M{r)]q{r) + q'{r)] 



dq {r) 1 .,, , .,, , 

-1-^ = q' {r) + lup' {r) 

dr r 



(A8) 



x{[e'-M(r)]p(r) +p'{r)} 
+ ^^T^,[e'-M{r)]p{r). 

In the nonrelativistic limit (c — >■ 00), these formulas re- 
duce to Eqs. ([HI), ([HI), (HID, and 
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